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Abstract 

We present an analytical study of the time dependent diffusion coefficient in a dilute suspension 
of spheres with partially absorbing boundary condition. Following Kirkpatrick (J. Chem. Phys. 
76, 4255) we obtain a perturbative expansion for the time dependent particle density using volume 
fraction / of spheres as an expansion parameter. The exact single particle i-operator for partially 
absorbing boundary condition is used to obtain a closed form time-dependent diffusion coefficient 
D{t) accurate to first order in the volume fraction /. Short and long time limits of D{t) are checked 
against the known short-time results for partially or fully absorbing boundary conditions and long- 
time results for reflecting boundary conditions. For fully absorbing boundary condition the long 
time diffusion coefficient is found to be D(t) = ha 2 / {I2f D$t) + 0((Dot/a 2 )~ 2 ), to the first order 
of perturbation theory. Here / is small but non-zero, Dq the diffusion coefficient in the absence of 
spheres, and a the radius of the spheres. The validity of this perturbative result is discussed. 
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I. INTRODUCTION 



There are numerous processes in biology, physics, geophysics, chemical engineering and 
materials science where a diffusing species reacts upon contact with the surface of another 
phase (typically a solid) and the effects of simultaneous diffusion and reaction are important. 
Studying the time-dependence of concentration of the species reacting with a porous host via 
diffusion, which goes back to Smoluchowsky (1913) x , remains an important problem both for 
practical and theoretical reasons: for reviews, see Weiss J , Kayser and Hubbard ', Fixman 4 , 
Torquato ', deGennes'*. Heterogeneous catalysis, transport and absorption of nutrients and 
drugs are well known examples. Another example is relaxation of Nuclear magnetization 
of interstitial fluid by interaction with walls in the NMR experiments. The connection 
between surface relaxation and NMR goes back to the seminal 1951 paper "On nuclear 
relaxation in gases by surface catalysis" by Bloch' himself. NMR Relaxation is now routinely 
used in numerous settings ranging from clinical MRF to geological explorations for energy 9 . 
The purpose of this paper is to investigate time- dependence of the diffusion coefficient of 
molecules confined in a well-connected porous host material with reactive walls. 

Diffusion measurements on surviving number of random walkers can give deep insight 
into this complex process of diffusion-relaxation. Experimentally observing the dynamics 
(diffusion) of molecules reveals a distinct set of information that cannot be obtained from 
the traditional relaxation (i.e. number of surviving walkers), measurements 11 . The complex 
porous systems, in whose interstices the fluid reside, have in general, convoluted and exten- 
sive well connected pores. The time dependence of the rate of relaxation in a well connected 
geometry which does not have arbitrarily large voids (see below) is not fully understood, 
not to mention diffusion. The relaxation processes are well understood for isolated cells or 
pores (isolated cavities), where an eigen-mode decomposition 12 ' 13 ' 14 of the diffusion equation 
with partially absorbing boundary condition is possible. Disentangling the complexity due 
to the combined effects of disorder, diffusion and decay remains a major hurdle in connected 
porous media 11 . 

A general theory of time-dependent diffusion and relaxation in porous media does not 
exist. At short-times, there are asymptotic results that are robust for all boundary conditions 
and arbitrary geometry 11 ' 15 ' 16 . For long-time, however, no general technique exists, except 
for simple models of a dilute periodic array of spheres with weak absorption 1 ' . In random 



2 



systems, a dilute suspension of randomly placed spheres is most well studied. For reflecting 
(non-absorbing) boundary condition de Swiet and Sen ls worked out a perturbation expansion 
based on Bixon and Zwanzig 1 " that gives short and long time behavior of the time dependent 
diffusion coefficient. 

The main goal of this paper is to extend these results 18 for partially absorbing boundary 
conditions. For partially absorbing boundary conditions perturbation expansions fail at 
long times 10 ' 1 '. The long time behavior of diffusion in such a system has not been previously 
addressed either analytically or numerically. Decay in particle number density with fully 
absorbing boundary condition has been studied extensively, see for example the review by 
Weiss 2 . The simplest effective medium approach of Smoluchowsky gives an exponential 
decay for the total number of surviving particles N(t). Bixon and Zwanzig pioneered the 
perturbative approach which gave a surprisingly slow power law decay. This discrepancy 
has triggered an intense interest in the problem. Kirkpatrick 10 perfected the perturbation 
theory, and he showed that by summing a select group of terms (the most divergent) this 
power law tail was removed and the exponential decay was recovered. However, the higher 
order terms in the perturbation series diverge as well and are difficult to re-sum and one 
cannot rule out long-time power law tails'". 

Later, allowing arbitrarily large fluctuations in local density, Grassberger and Proccacia, 
Kayser and Hubbard, and others 3,20 ' 21 conclusively showed that the long time limiting be- 
havior of N(t) in such system shows a non-perturbative behavior: in a <i-dimension system 
N(t) shows a stretched exponential decay N(t) — > exp[— £ d /( d + 2 )] as t — > oo. 21 This result 
cannot be captured by any (known) re-summation of the perturbation expansion. However 
three issues are important — (i) real systems do not have large voids and this stretched expo- 
nential result is of limited use (ii) Fixman 1 has argued that the crossover from the effective 
medium theory (simple exponential decay in time) to the exact limiting behavior (i.e. the 
stretched exponential) can be extraordinarily slow, i.e. the breakdown occurs at so long a 
time that N(t) is too small to be of any experimental relevance, (iii) we are interested in 
D(t) and not, per se, in N(t). D(t) is given by mean-square displacement divided by N(t), 
there is no known rigorous bound for this ratio. It is possible that perturbative expansion 
can give a reasonable answer for D(t). We note in passing that Grassberger and Proccacia 20 
give a heuristic argument that D(t) approaches zero as £- d /( d + 2 ) at long times. Even such a 
heuristic conjecture does not exist for partially absorbing systems nor for systems that do 
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not have large fluctuations in density. 

Therefore, we may argue that, for experimentally accessible times t, the perturbation 
expansion describes diffusion reasonably well and this is what we proceed to carry out. In 
the next section II we carry out a perturbation expansion for a dilute suspension of randomly 
placed spheres. We obtain an analytical expression of the approximate Green's function from 
first order terms of that expansion. From the Green's function we find the simple short and 
long time limiting behavior of time dependent diffusion coefficient D(t). 



II. ANALYTICAL RESULTS FOR SUSPENSION OF SPHERE 

In this section we study the problem of diffusion in dilute system with the partially 
absorbing boundary condition analytically. Only dilute systems are amenable to analytical 
methods. We will address the range of validity of such methods at the end of this section. 

Let us consider now diffusion of particles around a randomly distributed dilute suspension 
of N spheres, all of which of a given radius a, centered at positions Ri, 1 < % < N. The 
density C(f,t) of particles is governed by diffusion equation 

^ = D^C(r,t) (1) 

and the partially absorbing boundary condition is imposed: 

Dom-VCdf-Ril =a,t)=pC(\f-Ri\ =a,t), (2) 

C{\r-Ri\<a,t) = 0, (3) 

for % from 1 to N, where Hi is the unit vector normal to the ith sphere, pointing outward. 
The initial condition is 

C(f,t = 0) = C (f). ( 4 ) 
Following Kirkpatrick 1 " we go into the Laplace domain for time and the Fourier domain for 
space and express density in terms of i-operator of individual sphere as scattering center, 

N 

C(q,e) = G (q)C (q) + Y,G (q)T i ($C i (q,e), (5) 

i=i 

N 

G&e) = G (q)C (q) + ^ G^TffiCfa e). (6) 
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Here e is the reciprocal Laplace variable for time and q is the reciprocal Fourier variable 
for space and Go(q) = (e + Dq 2 )' 1 is the free particle Green's function, i.e. diffusion in 
unbounded space, and Tj is the Fourier space t-operator of a single sphere centered at 
Ri. Combining Eq. (5) and Eq. (6) we have the following perturbative expansion for the 
density 10 : 



{N N N 

G (q) +J2G (q)Ti(q)G (q) + J2J2 G ^ G ^ T MGo(q)+ 

N N N 1 

E E E GvT^G^T^G^TMGM + •■■ >C (<f), 

i=l j^i k^j ) 

This equation has the same form of the binary collision expansion in kinetic theory of gases 10 . 
This naive expansion encounters a problem at long time 10 ' 22 and it is better to sum the series 
to avoid such difficulties 10,19 . The averaged results reads: 

C(q,e) = {e + /V-5> 4 / J]dil t B^q, R x . . . R\, e)} C (q). (8) 

i=i J j=i 

Here n = N/V is the number of spheres per unit volume and Bi(q, Ri, e) represents the "self 
energy" operator for multiple scattering and includes correlation effects of averaging of i 
scattering center. For example, when Eq. (8)is expanded, the operator form gives 

J dRxB^q, Ri, e) = J dR^q), (9) 



dR 1 dR 2 B 2 (q,R 1 ,R 2 ,e)= / dR\ dR 2 {T 1 (q)Go(q)T 2 (q)G (q)T 1 (q) + 

J J (10) 

TMGo{q)TMG {q)TMGo{q)T 2 {q) + •■•}. 

As usual, here Tj is the t-operator of a single sphere at position Ri and the integrals are 
averaging over trap positions. In this paper we confine ourself to the effect of the first term 
of the expansion, i.e. we include multiple scattering from a single sphere to all orders, but 
do not consider correlation effects of multiple spheres. 

As required by translational invariance, the dependence of the t-operator on the position 
of sphere Ri is very simple: 

ti(q, q>, e, Ri) = exp[-iRi ■ (q- q')]t(q, qf,e). (11) 



Therefore the lowest order position averaged t-operator simply gives a factor of (27r) 3 5 3 (q— qf) 
and we only need the diagonal element of the t-operator (t) = t(q, q, e). 
Transformed into Laplace space the diffusion equation reads: 



(D V 2 -e)C(r,e) = -C (r), (12) 
The boundary condition in Eq. (2) becomes, in Laplace domain: 

Xani ■ VC(f, e) - C(f, e) = 0; A = — (13) 

pa 

on the surface \f — Ri\ = a. Here we introduced a dimensionless parameter A = Do/pa 
as a measure of strength of absorption on surface. A = corresponds to totally absorbing 
boundary condition and A = oo corresponds to reflecting boundary condition. 

It is convenient to choose the initial condition Co(r) = 5 3 (r — ff), which when Fourier 
transformed against ri) reads: 

e i<f/-f 

^o(r,<M = — . (14) 

A special solution for inhomogeneous equation (12) is 

1 e^'-r 

tf^(^ C )= (2w) 3 e + Dog/a - ( 15 ) 

and the general solution of the homogeneous equation satisfying boundary condition that 
density to be finite at infinity is: 

Cgeneral(r,e) = ybjP^COsO')^ ), (16) 

i=o a 

where Ki(x) is the ith order modified spherical Bessel function and Pi(x) is the zth order 
Legendre polynomial. 9' is the angle between <f and r — Ri and we introduce dimensionless 
parameter s = ea 2 /D to replace e. 

Choosing constants &j of Eq. (16) to make a linear combination with Eq. (15) that satisfy 
boundary condition Eq. (13) we have for \f — Ri\ > a: 



C(f, qf, e) 



(2tt) 3 e + D qf 2 



Jq/-(r-Ri) 



oo 

■E 

n=0 



{2n+l)i n P n (cos9') 



i.Kq'^j'niq' ) -Jn(q'a) ,\r-Ri\y/s 



(17) 



6 



where j n (x) is the nth order spherical Bessel function. 

Fourier transforming Eq. (17) with respect to r we can easily obtain the single particle t 
operator: 



t(q, qf, e) 



\k>f n {k')- 3n (k>) 



iTraDoe^-^-^ ) s + k 2 



(27T) 



\k-k/\ 



J 1 {\k-k/\) + J2( 2n + ^)Pn{cOs6 



IX 



n=0 



(18) 



[Vs7n(^)«n+i(V«) - kj n+1 (k)K n (yG)] \; k = qa; kl = qla 



Here we made the obvious choice for dimensionless momentum k = qa and kf = qla and 9 is 
the angle between k and kf. This form for the t-operator manifestly obeys Eq. (11). Setting 
k = kf we easily obtain the diagonal element of the t-operator 



(t) = -(47raD )\]- i (s + k 2 ) 



?1=0 



n+1) 



W n (k)-j n (k) 



[\fsjn{k)n n+ i(^fs) - kj n+1 (k)n n (y/s)] 

(19) 



To obtain the effective diffusion coefficient, we expand Eq. (19) to order k 2 and obtain 
the following leading order behavior: 



-k 



-1 + (3 + V5)A) + (1 + 2v/i + s)\ 2 



(t) = - (4 ^ 0) ^ (S+ " 2) + rTOTi " L3(l + A + A^)(l + v ^ + (2 + 2^ +S )A 



(20) 



Taking the lowest order in Eq. (8), combining with Eq. (9) and definition of Green's function 

Y-l N 

r ~~ V 



we have (G) = (G 1 — £{t}) 1 . Substituting from Eq. (20) we have, again to the order of 



k 2 



! + /(!" II 



-l+(3+ v /i)A+(l+ v /i) 2 A 2 



(l+A+A v ^)(l+ v / i+A(2+2 v / i+s)) < 



s + f(s + 



3(l+y^) ■ 
(1+A+A^) - 



Here we define the volume fraction of the spheres / = A ™y N . This generalizes Kirkpatrick's 1 
Eq. (A5) to finite absorption. To obtain D(t) notice that for the system with partially 



(21) 



r'olO 
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absorbing boundary condition: 

lim dG{k,t)/dk 2 
< r2 > = - 6a ""l imG(M ) ' (22) 

The numerator is the inverse Laplace transform of the following expression: 

-l+(3+V5)A+(l+y^) 2 A 2 



+ (1+A+AV5' -+y/5+A(2+2y5+s)) 

The denominator, that is the number of surviving particles N(t), is the inverse Laplace 
transform of G(k = 0): 

(24) 



Performing the inverse Laplace transform of the above equation, we obtain for the de 
nominator, the total number of surviving particles: 



N(t) 



(1 + f)(n - r 2 )(n - r 3 )(r 2 - r 3 )A 



e r ^r 1 (r 2 — r 3 )(r 1 A + A + l)erfc ^— T\\fzj + c.p. 

(25) 

here t = Dot /a 2 , erfc(x) = 1 — erf(x) is the complementary error function, c.p. denotes 
cyclical permutation n — > r 2 , r 2 — >• r 3 and r 3 — > ri, and ri,r 2 ,r 3 are the roots of the cubic 
equation: 

(1 + A)/r 3 + (1 + / + Xf)r 2 + 3/r + 3/ = 0. (26) 

To check this with the known results we look at several limits where the forms of N(t) are 
simple. 

First consider the short time limit t < 1, i, e., t a 2 / Dq. There are two subcases. If 



the boundary is nearly reflecting, specifically s/s 3> 1/A, i.e. p <^ y/D /t, the results is 
particularly simple: 

N(t) = l-f-^ + ^ + (t 5 / 2 ) . (27) 

Here 1 — / factor comes from the excluded volume of the spheres. The first term agrees with 
equation (A9) of Mitra et. al. w , which is only to the first order in /. 



In the other extreme, when p ^> \/D /t, that is the nearly fully absorbing boundary 
condition, the short time total number of particle becomes insensitive to p: 

N(t) = 1 - / - (6/ ~y )v/? + (-3/ + 15f 2 )t + 0(F/ 2 ), (28) 

'71 



The first term in turns agrees with equation (A8) of Mitra et. al. 16 , again to the first order 
in /. 

Next consider the long-time limits. The total particle number decays, in the long time 
limit, t 3> 1, i.e., t ^> a 2 /D , for fully absorbing boundary condition: 

mi) - ^ **" + 2 ~ f ^ f)X + o<r<*) (29, 

The first term is twice the perturbative result in Bixon and Zwanzig 19 . Recall that the exact 
non-perturbative result of Grassberger and Procaccia 20 gives a stretched exponential decay. 
Also note the 1/f dependence, for small but finite /. 

Finally consider the diffusion coefficient with partially absorbing boundary conditions. 
Results from the inverse Laplace transform of the numerator in Eq. (23) is too complex to 
be reproduced fully here, so we summarize the result in a simpler form: 

Ss = E E l ) + E *); ( 30 ) 

i=l j=l i=4 

In this equation, ri_ 3 are the same three roots of the cubic equation Eq. (26) and = 
( — 1 — 2A ± y/1 — 4A 2 )/2A are the two roots of the quadratic equation. And the constants 
Aij, which depends on r\ through r 5 above as well as A and /, are the coefficients of partial 
fractions of Eq. (23): 

[-I , f(-i _ -l+(3+y^)A+(l+Vi) 2 A 2 n-I 

L ~ v (l+A+A v ^)(l+ v ^+A(2+2 v ^+s))^J , s 

3(1+ Vi) m2 W 



[' + /('+75Sgfe)] 



5 3 A i2 

= S (v^-V- (32) 

7^2 are given by the following expression obtained by inverse Laplace transform of the partial 
fractions: 



7i(ri, t) = — = + Tie Ti erfc(-rjV £ 
Vvrt 

I 2 (r u t) = + (1 + 2r 2 )r i e^erfc(-r l Vt) (33) 

At short time t <C a 2 /D and for the reflecting boundary condition, using the fact N(t) = 
1 — /, we recover the results of Eq. (14)of De Swiet and Sen 18 : 

D(i)=D [l-¥jl + O®]. (34) 
9 



For the fully absorbing boundary condition at short time, we expand Eq. (30) to obtain: 

This combined with N(t) obtained in Eq. (28), gives to the first order of / 

D(t)/D = l-^ + O(P), (36) 

which is consistent with the result in Mitra et. al. 16 Eq. (All) 25 . 

Much more interesting is the case of fully absorbing boundary condition at long time 
t > a 2 /D . There the numerator of Eq. (23) reduced to {f(f(6f + 29) + 19) + 5}{18/ 2 (/ + 
1)(/ + 4) v /7r} t -3 / 2 + 0(t~ 5 / 2 ). Keeping / to first order and combining with Eq. (22) and 
the result for N(t) Eq. (29) we obtain the surprising result: 

D(t)/D = -^r 1 + O(r% (37) 

that is, within perturbation theory, the long time diffusion coefficient D(t) of the dilute 
suspension of spheres with absorbing boundary condition approaches zero as \jt. Notice 
that D(t) is inversely dependent on volume fraction / as is N(t). There is a subtle order of 
limit: though we are taking an asymptotic expansion for small /, the limit where Eq. (37) 
is valid is for a small yet fixed / but t — > oo, i.e. t > 1//, where the limit is well defined in 
the equation. 

Next consider some known results. For diffusion in a straight tube with fully absorbing 
wall, the separation of variables can be used to obtain exact results. The number of particles 
N(t) will decay exponentially to zero with time, but D(t) — > -Do/3 as t — > oo because the 
particles may diffuse along the axis of the tube. Any long tube-like open pore with slowly 
changing diameter will also have finite diffusion coefficient. Furthermore a periodic array of 
spheres with fully absorbing boundary condition always have finite diffusion coefficient at 
long time at any volume fraction. In fact there is no known connected system, excluding 
isolated-pore like structures (e.g. in the Lifshitz limit 20 where the absorbers form a cavity), 
with any amount of surface relaxation, which has a vanishing long time diffusion coefficient. 
Intuitively a dilute random suspension of spheres is not likely to be an exception. So the 
perturbative result of 1/t behavior of D(t) may indicate a rapid drop of diffusion coefficient 
at a certain range of t but is unlikely to give the correct true long time behavior. 
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III. CONCLUSION 



The 1/t time dependence for D(t) given by the perturbative result for random suspension 
is peculiar, as is the 1/t 3 / 5 behavior given by the heuristic arguments of Grassberger and 
Procaccia 20 . While analytical perturbative method considered here give deep insight and 
give correct result for short-time limit, we suspect that perturbative results for the long 
time behavior of D{t) will break down just as the perturbative result fails for total number 
of surviving particles N(t) at long times. In view of the findings of Fixman 1 , it will be 
interesting to examine the perturbative results for D(t) against numerical results with a 
hope that for the time-regime of experimental interest, perturbative results will suffice. 
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